use "C:\PATH\final_data_replication", clear
drop if t>36 | t<5




***Part 1: Their aggregation method is not robust to (a)correcting for the errors and/or (b) using the 2019 reform which is the only one with subsidy changes for both types of communities


cap: egen tmonopoly=group(t monopoly)


***Their specification
reghdfe rnfb rnfb_sub if t>4 & t<=36, absorb(community t) vce(cluster community)
gen sample=1 if e(sample)==1

reghdfe rnfb rnfb_sub if t>4 & t<=36, absorb(community tmonopoly) vce(cluster community)
gen qsubsidypost=qsubsidy if t>32

egen totalsubsidy=total(qsubsidy*e(sample))
egen totalsubsidypost=total(qsubsidypost*e(sample))
bysort monopoly: egen share=total((qsubsidy/totalsubsidy)*e(sample))
bysort monopoly: egen sharepost=total((qsubsidypost/totalsubsidypost)*e(sample))

bysort monopoly:  summ share sharepost

local sharenon=0.864
local sharemonop=0.136


***Full sampe original data (replicating Table 1 in their paper)

reghdfe rnfb rnfb_sub if t>4 & t<=36, absorb(community t) vce(cluster community)
outreg2 using regtable1a, tex label ctitle(rnfb) bdec(3) rdec(3) alpha(0.01,0.05,0.10) addstat(Adj R-squared,e(r2_a))	nocons noaster replace

reghdfe rnfb rnfb_sub rnfb_sub_monopoly if t>4 & t<=36, absorb(community t) vce(cluster community)
lincom (_b[rnfb_sub]*`sharenon')+((_b[rnfb_sub]+_b[rnfb_sub_monopoly])*`sharemonop')
local j2=r(ub)
local j1=r(lb)
local j0=r(estimate)
local fmtja: display  %4.3f `j0'
local fmtj2a: display %4.3f `j2'
local fmtj1a: display %4.3f `j1'
outreg2 using regtable1a, tex label ctitle(rnfb) bdec(3) rdec(3) alpha(0.01,0.05,0.10) addstat(Adj R-squared,e(r2_a))	addtext(Weighted Avg., `fmtja', 95% CI,"(`fmtj1a', `fmtj2a')") nocons noaster append

reghdfe rnfb rnfb_sub rnfb_sub_monopoly if t>4 & t<=36, absorb(community tmonopoly) vce(cluster community)
lincom (_b[rnfb_sub]*`sharenon')+((_b[rnfb_sub]+_b[rnfb_sub_monopoly])*`sharemonop')
local j2=r(ub)
local j1=r(lb)
local j0=r(estimate)
local fmtja: display  %4.3f `j0'
local fmtj2a: display %4.3f `j2'
local fmtj1a: display %4.3f `j1'
outreg2 using regtable1a, tex label ctitle(rnfb) bdec(3) rdec(3) alpha(0.01,0.05,0.10) addstat(Adj R-squared,e(r2_a))	addtext(Weighted Avg., `fmtja', 95% CI,"(`fmtj1a', `fmtj2a')") nocons noaster append


*twoway (scatter coeffm year `restrict', yline(0) ytitle("Coefficient") title("Hunt occasionally/frequently: Food Mail-NNC") lcolor(black) legend(off) xline(2011)) (rcap upcifm locifm year))

reghdfe rnfb rnfb_sub if t>=25 & t<=36, absorb(community t) vce(cluster community)
outreg2 using regtable1a, tex label ctitle(rnfb) bdec(3) rdec(3) alpha(0.01,0.05,0.10) addstat(Adj R-squared,e(r2_a))	nocons noaster append
reghdfe rnfb rnfb_sub rnfb_sub_monopoly if t>=25 & t<=36, absorb(community t) vce(cluster community)
lincom (_b[rnfb_sub]*`sharenon')+((_b[rnfb_sub]+_b[rnfb_sub_monopoly])*`sharemonop')
local j2=r(ub)
local j1=r(lb)
local j0=r(estimate)
local fmtja: display  %4.3f `j0'
local fmtj2a: display %4.3f `j2'
local fmtj1a: display %4.3f `j1'
outreg2 using regtable1a, tex label ctitle(rnfb) bdec(3) rdec(3) alpha(0.01,0.05,0.10) addstat(Adj R-squared,e(r2_a))	addtext(Weighted Avg., `fmtja', 95% CI,"(`fmtj1a', `fmtj2a')") nocons noaster append
reghdfe rnfb rnfb_sub rnfb_sub_monopoly if t>=25 & t<=36, absorb(community tmonopoly) vce(cluster community)
lincom (_b[rnfb_sub]*`sharenon')+((_b[rnfb_sub]+_b[rnfb_sub_monopoly])*`sharemonop')
local j2=r(ub)
local j1=r(lb)
local j0=r(estimate)
local fmtja: display  %4.3f `j0'
local fmtj2a: display %4.3f `j2'
local fmtj1a: display %4.3f `j1'
outreg2 using regtable1a, tex label ctitle(rnfb) bdec(3) rdec(3) alpha(0.01,0.05,0.10) addstat(Adj R-squared,e(r2_a))	addtext(Weighted Avg., `fmtja', 95% CI,"(`fmtj1a', `fmtj2a')") nocons noaster append

reghdfe rnfb rnfb_sub if t>=29 & t<=36, absorb(community t) vce(cluster community)
outreg2 using regtable1a, tex label ctitle(rnfb) bdec(3) rdec(3) alpha(0.01,0.05,0.10) addstat(Adj R-squared,e(r2_a))	nocons noaster append
reghdfe rnfb rnfb_sub rnfb_sub_monopoly if t>=29 & t<=36, absorb(community t) vce(cluster community)
lincom (_b[rnfb_sub]*`sharenon')+((_b[rnfb_sub]+_b[rnfb_sub_monopoly])*`sharemonop')
local j2=r(ub)
local j1=r(lb)
local j0=r(estimate)
local fmtja: display  %4.3f `j0'
local fmtj2a: display %4.3f `j2'
local fmtj1a: display %4.3f `j1'
outreg2 using regtable1a, tex label ctitle(rnfb) bdec(3) rdec(3) alpha(0.01,0.05,0.10) addstat(Adj R-squared,e(r2_a))	addtext(Weighted Avg., `fmtja', 95% CI,"(`fmtj1a', `fmtj2a')") nocons noaster append
reghdfe rnfb rnfb_sub rnfb_sub_monopoly if t>=29 & t<=36, absorb(community tmonopoly) vce(cluster community)
lincom (_b[rnfb_sub]*`sharenon')+((_b[rnfb_sub]+_b[rnfb_sub_monopoly])*`sharemonop')
local j2=r(ub)
local j1=r(lb)
local j0=r(estimate)
local fmtja: display  %4.3f `j0'
local fmtj2a: display %4.3f `j2'
local fmtj1a: display %4.3f `j1'
outreg2 using regtable1a, tex label ctitle(rnfb) bdec(3) rdec(3) alpha(0.01,0.05,0.10) addstat(Adj R-squared,e(r2_a))	addtext(Weighted Avg., `fmtja', 95% CI,"(`fmtj1a', `fmtj2a')") nocons noaster append




*Fixing errors

**Error 1: mis-classification of Clyde River as non-monopoly
replace monopoly=1 if community=="Clyde River (Kangiqtugaapik)"


cap: drop tmonopoly
cap: egen tmonopoly=group(t monopoly)

replace rnfb_sub_monopoly=rnfb_sub*monopoly

**need to recalculate the shares
bysort monopoly: egen sharefix=total((qsubsidy/totalsubsidy)*sample)
bysort monopoly: summ sharefix
local sharenon=0.838
local sharemonop=0.162

**Error 2: Stony Rapids and Black Lake did not ever receive the subsidy
*Note that these communities do not have any price observations after t=30
gen problem=0
replace problem=1 if community=="Stony Rapids" | community=="Black Lake"
replace rnfb_sub=0 if problem==1
replace rnfb_sub_monopoly=0 if problem==1



**Full sample, corrected data
reghdfe rnfb rnfb_sub if t>4 & t<=36, absorb(community t) vce(cluster community)
outreg2 using regtable1b, tex label ctitle(rnfb) bdec(3) rdec(3) alpha(0.01,0.05,0.10) addstat(Adj R-squared,e(r2_a))	nocons noaster replace

reghdfe rnfb rnfb_sub rnfb_sub_monopoly if t>4 & t<=36, absorb(community t) vce(cluster community)
lincom (_b[rnfb_sub]*`sharenon')+((_b[rnfb_sub]+_b[rnfb_sub_monopoly])*`sharemonop')
local j2=r(ub)
local j1=r(lb)
local j0=r(estimate)
local fmtja: display  %4.3f `j0'
local fmtj2a: display %4.3f `j2'
local fmtj1a: display %4.3f `j1'
outreg2 using regtable1b, tex label ctitle(rnfb) bdec(3) rdec(3) alpha(0.01,0.05,0.10) addstat(Adj R-squared,e(r2_a))	addtext(Weighted Avg., `fmtja', 95% CI,"(`fmtj1a', `fmtj2a')") nocons noaster append

reghdfe rnfb rnfb_sub rnfb_sub_monopoly if t>4 & t<=36, absorb(community tmonopoly) vce(cluster community)
lincom (_b[rnfb_sub]*`sharenon')+((_b[rnfb_sub]+_b[rnfb_sub_monopoly])*`sharemonop')
local j2=r(ub)
local j1=r(lb)
local j0=r(estimate)
local fmtja: display  %4.3f `j0'
local fmtj2a: display %4.3f `j2'
local fmtj1a: display %4.3f `j1'
outreg2 using regtable1b, tex label ctitle(rnfb) bdec(3) rdec(3) alpha(0.01,0.05,0.10) addstat(Adj R-squared,e(r2_a))	addtext(Weighted Avg., `fmtja', 95% CI,"(`fmtj1a', `fmtj2a')") nocons noaster append





***Use only 2019 data (next 3)--> same as in the original paper...

reghdfe rnfb rnfb_sub if t>=25 & t<=36, absorb(community t) vce(cluster community)
outreg2 using regtable1b, tex label ctitle(rnfb) bdec(3) rdec(3) alpha(0.01,0.05,0.10) addstat(Adj R-squared,e(r2_a))	nocons noaster append

reghdfe rnfb rnfb_sub rnfb_sub_monopoly if t>=25 & t<=36, absorb(community t) vce(cluster community)
lincom (_b[rnfb_sub]*`sharenon')+((_b[rnfb_sub]+_b[rnfb_sub_monopoly])*`sharemonop')
local j2=r(ub)
local j1=r(lb)
local j0=r(estimate)
local fmtja: display  %4.3f `j0'
local fmtj2a: display %4.3f `j2'
local fmtj1a: display %4.3f `j1'
outreg2 using regtable1b, tex label ctitle(rnfb) bdec(3) rdec(3) alpha(0.01,0.05,0.10) addstat(Adj R-squared,e(r2_a))	addtext(Weighted Avg., `fmtja', 95% CI,"(`fmtj1a', `fmtj2a')") nocons noaster append

reghdfe rnfb rnfb_sub rnfb_sub_monopoly if t>=25 & t<=36, absorb(community tmonopoly) vce(cluster community)
lincom (_b[rnfb_sub]*`sharenon')+((_b[rnfb_sub]+_b[rnfb_sub_monopoly])*`sharemonop')
local j2=r(ub)
local j1=r(lb)
local j0=r(estimate)
local fmtja: display  %4.3f `j0'
local fmtj2a: display %4.3f `j2'
local fmtj1a: display %4.3f `j1'
outreg2 using regtable1b, tex label ctitle(rnfb) bdec(3) rdec(3) alpha(0.01,0.05,0.10) addstat(Adj R-squared,e(r2_a))	addtext(Weighted Avg., `fmtja', 95% CI,"(`fmtj1a', `fmtj2a')") nocons noaster append


***********
reghdfe rnfb rnfb_sub if t>=29 & t<=36, absorb(community t) vce(cluster community)
outreg2 using regtable1b, tex label ctitle(rnfb) bdec(3) rdec(3) alpha(0.01,0.05,0.10) addstat(Adj R-squared,e(r2_a))	nocons noaster append

reghdfe rnfb rnfb_sub rnfb_sub_monopoly if t>=29 & t<=36, absorb(community t) vce(cluster community)
lincom (_b[rnfb_sub]*`sharenon')+((_b[rnfb_sub]+_b[rnfb_sub_monopoly])*`sharemonop')
local j2=r(ub)
local j1=r(lb)
local j0=r(estimate)
local fmtja: display  %4.3f `j0'
local fmtj2a: display %4.3f `j2'
local fmtj1a: display %4.3f `j1'
outreg2 using regtable1b, tex label ctitle(rnfb) bdec(3) rdec(3) alpha(0.01,0.05,0.10) addstat(Adj R-squared,e(r2_a))	addtext(Weighted Avg., `fmtja', 95% CI,"(`fmtj1a', `fmtj2a')") nocons noaster append

reghdfe rnfb rnfb_sub rnfb_sub_monopoly if t>=29 & t<=36, absorb(community tmonopoly) vce(cluster community)
lincom (_b[rnfb_sub]*`sharenon')+((_b[rnfb_sub]+_b[rnfb_sub_monopoly])*`sharemonop')
local j2=r(ub)
local j1=r(lb)
local j0=r(estimate)
local fmtja: display  %4.3f `j0'
local fmtj2a: display %4.3f `j2'
local fmtj1a: display %4.3f `j1'
outreg2 using regtable1b, tex label ctitle(rnfb) bdec(3) rdec(3) alpha(0.01,0.05,0.10) addstat(Adj R-squared,e(r2_a))	addtext(Weighted Avg., `fmtja', 95% CI,"(`fmtj1a', `fmtj2a')") nocons noaster append





*****************************Part 2*******************
**Show that the nature of the subsidy variation used matters (within/across monopolies) including pre/post estimation
**Discuss issues of statistical power and multiple testing in the context of heterogeneous effects


bysort community year: egen yrnfb=mean(rnfb)

gen post=1 if t>=33
gen pre=1 if t>=25 & t<33
bysort community: egen postrnfb=mean(rnfb*post)
bysort community: egen prernfb=mean(rnfb*pre)
gen dpostpre=postrnfb-prernfb

cap: drop dprice dsubsidy
sort community t
by community: gen dprice4=rnfb-rnfb[_n-4] if t==t[_n-4]+4
by community: gen dyprice4=yrnfb-yrnfb[_n-4] if t==t[_n-4]+4

replace dr2019=0 if problem==1

gen mdr2019=-dr2019



***adjust axes
***this does not drop any of the observed price changes but makes the axes nice
drop if dr2019>34
drop if dr2019<12
 
summ dr2019 if rnfb~=. & dr2019~=0 & year==2019, detail

set scheme s1color
**Using Q4 only
twoway (scatter dprice4 dr2019 if t==36 & monopoly==1, ytitle("Change in price") xtitle("Change in subsidy") title(2019Q4 vs. 2018Q4) mcolor(red) legend(order(1 "Monopoly" 2 "Non-Monop." 3 "Monopoly" 4 "Non-monop." 5 "Pooled" 6 "dP=dS"))) (scatter dprice4 dr2019 if t==36 & monopoly==0, mcolor(black))  (lfit dprice4 dr2019 if t==36 & monopoly==1, lcolor(red))  (lfit dprice4 dr2019 if t==36 & monopoly==0, lcolor(black)) (lfit dprice4 dr2019 if t==36, lcolor(blue)) (lfit mdr2019 dr2019 if t==36 & dr2019>16 & dr2019<34, lcolor(green)) 
graph save leftpanel, replace


**Using annual averages for 2018 and 2019
twoway (scatter dyprice4 dr2019 if t==36 & monopoly==1,  ytitle("Change in price") xtitle("Change in subsidy") title(2019 vs. 2018) mcolor(red) legend(order(1 "Monopoly" 2 "Non-monop." 3 "Monopoly" 4 "Non-monop." 5 "Pooled" 6 "dP=dS"))) (scatter dyprice4 dr2019 if t==36 & monopoly==0, mcolor(black))  (lfit dyprice4 dr2019 if t==36 & monopoly==1, lcolor(red))  (lfit dyprice4 dr2019 if t==36 & monopoly==0, lcolor(black)) (lfit dyprice4 dr2019 if t==36, lcolor(blue)) (lfit mdr2019 dr2019 if t==36, lcolor(green)) 
graph save rightpanel, replace


set scheme s2color
graph combine leftpanel.gph rightpanel.gph, ycommon xcommon
graph export figure1.png, replace
graph export figure1.pdf, replace


***These graphs make the point about the nature of the variation?




***Power analysis
cap: drop dr
gen dr=dr2016 if t==24
replace dr=dr2019 if t==36



*ex-ante power analysis: full sample
summ dr2019 dyprice4 if t==36 & dyprice4~=. & dr2019~=., detail
power oneslope -1, n(78) power(0.8) sdx(5.10) sdy(7.42) onesided
power oneslope -1 -0.66, n(78) sdx(5.10) sdy(7.42) onesided
power oneslope -1 -0.5, n(78) sdx(5.10) sdy(7.42) onesided


*Non-monopoly
summ dr2019 dyprice4 if t==36 & dyprice4~=. & dr2019~=. & monopoly==0, detail
power oneslope -1, n(45) power(0.8) sdx(4.28) sdy(7.12) onesided
power oneslope -1 -0.66, n(45) sdx(4.28) sdy(7.12) onesided
power oneslope -1 -0.5, n(45) sdx(4.28) sdy(7.12) onesided


*Monopoly
summ dr2019 dyprice4 if t==36 & dyprice4~=. & dr2019~=. & monopoly==1, detail
power oneslope -1, n(33) power(0.8) sdx(5.58) sdy(6.09) onesided
power oneslope -1 -0.66, n(33) sdx(5.58) sdy(6.09) onesided
power oneslope -1 -0.5, n(33) sdx(5.58) sdy(6.09) onesided










****Figure 1********
*copy-paste from spreadsheet
cap: drop counter
gen counter=_n
replace counter=counter-9 if counter>9


gen rcounter=10-counter if counter<10


label define counter 1 "(1) All periods pooled" 2 "(2) All periods NDK" 3 "(3) All periods NDK monop. x time FE" 4 "(4) 2017Q1-2019Q4 pooled" 5 "(5) 2017Q1-2019Q4 NDK" 6 "(6) 2017Q-2019Q4 NDK monop. x time FE" 7 "(7) 2018Q1-2019Q1 pooled" 8 "(8) 2018Q1-2019Q4 NDK" 9 "(9) 2018Q1-2019Q4 NDK monop. x time FE", replace

label values counter counter


label define rcounter 9 "(1) All periods pooled" 8 "(2) All periods NDK" 7 "(3) All periods NDK monop. x time FE" 6 "(4) 2017Q1-2019Q4 pooled" 5 "(5) 2017Q1-2019Q4 NDK" 4 "(6) 2017Q-2019Q4 NDK monop. x time FE" 3 "(7) 2018Q1-2019Q1 pooled" 2 "(8) 2018Q1-2019Q4 NDK" 1 " (9) 2018Q1-2019Q4 NDK monop. x time FE", replace

label values rcounter rcounter

set scheme s1mono
local restrict if sample==1
* & counter<10




*twoway (scatter coefficient counter `restrict', xlabel(1(1)9) ytitle("Pass-through coefficient") xtitle("Table 1 Panel A Column") title("Original data set") lcolor(black) legend(off) yline(-1) yline(-0.67, lpattern(dash))) (rcap ub lb counter `restrict')




twoway (scatter rcounter coefficient `restrict', ylabel(1 2 3 4 5 6 7 8 9, valuelabel angle(horizontal)) xtitle("Pass-through coefficient") ytitle("Table 1 Panel A Column") title("Original data set") lcolor(black) legend(off) xline(-1) xline(-0.67, lpattern(dash))) (rcap ub lb rcounter `restrict', horizontal)

graph save ndk1.gph, replace





local restrict if sample==2
* & counter>9
*twoway (scatter coefficient counter `restrict', xlabel(1(1)9) ytitle("Pass-through coefficient") xtitle("Table 1 Panel B Column") title("Corrected data set") lcolor(black) legend(off) yline(-1) yline(-0.67,lpattern(dash))) (rcap ub lb counter `restrict')

twoway (scatter rcounter coefficient `restrict', ylabel(1 2 3 4 5 6 7 8 9, valuelabel angle(horizontal)) xtitle("Pass-through coefficient") ytitle("Table 1 Panel B Column") title("Corrected data set") lcolor(black) legend(off) xline(-1) xline(-0.67, lpattern(dash))) (rcap ub lb rcounter `restrict', horizontal)


graph save ndk2.gph, replace

graph combine ndk1.gph ndk2.gph, col(1) xcommon
graph export "C:\Users\nickl\Dropbox\nutrition north\comment on naylor\ndk.pdf", replace
graph export "C:\Users\nickl\Dropbox\nutrition north\comment on naylor\ndk.jpg", replace